Droplet evolution in expanding flow of warm dense matter 



Julien Armijcj^] 

Lawrence Berkeley National Laboratory, Berkeley, California, USA 

John J. Barnard 
Lawrence Livermore National Laboratory, Livermore, California, USA 
(Dated: January 26, 2013) 

We propose a simple, self-consistent kinetic model for the evolution of a mixture of droplets and 
vapor expanding adiabatically in vacuum after rapid, almost isochoric heating. We study the evo- 
lution of the two-phase fluid at intermediate times between the molecular and the hydrodynamic 
scales, focusing on out-of-equilibrium and surface effects. We use the van der Waals equation of state 
as a test bed to implement our model and study the phenomenology of the upcoming second neutral- 
ized drift compression experiment (NDCX-II) at Lawrence Berkeley National Laboratory (LBNL) 
that uses ion beams for target heating. We find an approximate expression for the temperature dif- 
ference between the droplets and the expanding gas and we check it with numerical calculations. 
The formula provides a useful criterion to distinguish the thermalized and nonthermalized regimes 
of expansion. In the thermalized case, the liquid fraction grows in a proportion that we estimate 
analytically, whereas, in case of too rapid expansion, a strict limit for the evaporation of droplets is 
derived. The range of experimental situations is discussed. 

PACS numbers: 64.70.fm, 51.10.+y 



I. INTRODUCTION 

Warm dense matter (WDM) conditions, corresponding 
roughly to densites 0.01 < p/paoUd < 10 and tempera- 
tures O.leV < T < 10 eV, can be defined as the region 
of thermodynamic space corresponding to the double 
crossover from degenerate to non-degenerate and from 
weakly to strongly coupled matter [T] , so that the "easy" 
limiting descriptions in terms of cool plasma and hot con- 
densed matter meet and have to be somehow connected 
to each other. This problem is drawing growing attention 
because of the serious theoretical challenges involved, and 
because of the occurence of WDM in the contexts of iner- 
tial fusion energy (IFE) , astrophysics (planet cores) , and 
laser ablation for materials processing, nanoparticles for- 
mation and film deposition pHl]- 

Generally, WDM experiments are inertially confined 
and explosive. Rapid heating of the material is re- 
quired, so that the energy deposition (by lasers, ions, 
neutrons, electrical discharges, etc.) is faster than its 
release through hydrodynamic expansion. Pressures in 
the kbar to Mbar range can be reached before giving rise 
to supersonic expansion with typical outflow velocities of 
several km/s. 

The two-phase region of the phase diagram belongs 
only partly to the WDM regime, with its high temper- 
ature part near the critical point. However, any WDM 
experiment will almost inevitably lead to two-phase con- 
ditions during the expansion and cooling process. This 
happens either from below the critical point (ion heating 
experiments at Gesellschaft fiir Schwerionenforschung in 
Darmstadt [SI |6], or second neutralized drift compres- 
sion experiment (NDCX-II) at Lawrence Berkeley Na- 
tional Laboratory (LBNL) f71, low fluence laser ablation. 



Z machines), or from above it (IFE, high fluence laser ab- 
lation, upcoming ion heating machines) . In the first case 
an overstretched liquid fragments and evaporates into a 
mixture of droplets and gas, whereas in the second case a 
hot supersaturated gas nucleates small clusters while ex- 
panding. In both cases the flow becomes a plume of gas 
and condensed clusters (most often liquid, so the term 
"droplet" is appropriate) : the monophasic liquid or gas 
has undergone phase separation with the creation of sur- 
faces giving a non-trivial geometry to the fluid, which 
may a priori affect its dynamical properties. 

Recently, there has been significant progress in the ob- 
servation of those two-phase flows, from the early abla- 
tion and plume expansion stages in the the ps and ns 
timescales [SUIT] to the late /is timescale evolution in- 
cluding postmortem analysis of the clusters [121 HI] • 

Basic questions arise when considering a two-phase 
flow. First, what is the droplets' size and distribution, 
and how do they evolve during the expansion? Second, 
are the droplets and the gas in thermal equilibrium? The 
answer can determine the conditions of validity for hydro- 
dynamic approaches based on the Maxwell construction 
or any two-phase equation of state (EoS) that assumes 
local equilibrium. 

So far, two main approaches have been used. On 
one hand, molecular dynamics (MD) codes ^14-,--20, com- 
pute the dynamics of each particle separately, and have 
given powerful insight into the processes of fragmenta- 
tion, phase explosion, and the different mechanisms for 
ablation, but they are inherently limited to only treat 
the early times (< Ins, |21j). and with a small number 
of particles (~ lO''). On the other hand, hydrodynamic 
codes [131 l22H25| can model experiments completely, but 
they deal with mesoscopic fluid cells, and often rely on 
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crude approximations concerning the molecular and ki- 
netic processes involved. Complex hydrodynamic codes 
including a treatment of the kinetics of phase change pro- 
cesses and surface effects in each cell are under develop- 
ment [131 US] , but providing a complete and reliable de- 
scription of a whole WDM experiment is still a challenge. 

Better understanding of two-phase flows should be 
helpful for the preparation of experiments, including the 
diagnostics, and for the interpretation of the data. In 
IFE especially, the problem of debris dynamics is a cru- 
cial issue due to their impact on the optics and other 
components of the target chambers 117| . 

In this paper we propose an alternative approach to 
study two-phase flows in the cool or late time WDM sit- 
uations. Our model was initially conceived to predict 
the phenomenology of the upcoming target heating ex- 
periments with the NDCX II machine at LBNL where 
an ion beam will almost isochorically heat a thin metal- 
lic foil to temperatures of about leV. However, the model 
should apply to any two-phase flow. 

The core of our model is a self-consistent set of kinetic 
rate equations for the particle and energy fluxes between 
a droplet and the surrounding gas in an expanding La- 
grangian cell. This set of equations can be applied to any 
two-phase EoS. The computing cell is considered as part 
of a larger hydrodynamic code, but in this paper we only 
consider one cell containing one droplet. We also neglect 
several features that could be added in the near future. 

To implement the kinetic equations and explore the 
patterns of two-phase expansion, we use the van der- 
Waals fluid model that makes it possible to build a com- 
plete set of thermodynamic functions. We thus demon- 
strate the ability of the kinetic model to simulate 
nonequilibrium two-phase flows in the wide range be- 
tween themolecular and hydrodynamic scales. In par- 
ticular, we use it to distinguish the different regimes of 
two-phase expansion: on one side, quasi- or fully ther- 
malized; on the other side, nonthemalized.We show that 
this distinction depends on the initial target dimensions 
and the initial temperature. We then study those regimes 
analytically and numerically. 

II. BACKGROUND 

A. Expanding two-phase flows. Supercritical and 
subcritical cases 

The model that we propose lies at a mesoscopic scale 
between the molecular and hydrodynamic scales, so we 
need some preliminar assumptions. Our computing cell 
is considered as an elementary piece of a larger hydro- 
dynamic code describing an expanding flow. The lin- 
ear strain rate 77 characterizes the expansion of the cell 
L = Lo{l + rit), where Lq is the initial cell size. We deflne 
the hydrodynamic time ihydro = ri~^. In the following. 



we assume rapid heating (thcating < ^hydro) so that en- 
ergy deposition in the material is almost isochoric. For 
simplicity, we assume instantaneous energy deposition. 

To get insight into the global flow, it is interesting to 
review some analytical and numerical results. The self 
similar rarefaction wave (SSRW) is the solution [55] de- 
scribing the one-dimensional (ID) expansion of a perfect 
gas (semi-infinite at z < 0) of adiabatic coeficient 7 af- 
ter instantaneous uniform heating at the initial temper- 
ature Tg. Denoting Cq the sound speed in the fluid at Tq, 
the SSRW solution describes an outward expanding front 
travelling at the outflow velocity vq = 2cg/ (7 — 1), which 
is 3co for a perfect monoatomic gas, while the inward rar- 
efaction wave propagates at eg [21]. Note that the SSRW 
can be computed semi-numerically for any EoS of a non- 
ideal gas [22] and has been validated by MD simulations 

m- 

As an example of a numerical simulation of expanding 
flows. Fig. [1] shows a hydrodynamic calculation with the 
code DPC using an EoS based on Maxwell construction 
[25] . Here the liquid and gas are assumed in equilibrium, 
which is not kinetically justified (see Sec. IV. B), and the 
outfiow velocity is about 8km/s after 10ns. 
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FIG. 1. Hydrodynamic calculation with DPC code of NDCX 
II reference case, from [24]. A 3.5/im-thick Al foil is heated 
within Ins with an ion beam and subsequently cools down 
during adiabatic expansion 

In the following we assume a flow with linear speed 
profile and outflow velocity Vq = 3cg, but it is worth re- 
marking that this is quite simplistic. In particular, sev- 
eral numerical works [22] [211 (25] 122] have reported the 
occurrence of "plateaus", that is, zones of nearly con- 
stant density, related to the fluid zones entering into the 
two-phase regime. 

Let us now present our simple classiflcation of two- 
phase expanding flows, which is based on equilibrium 
thermodynamics, and is similar to the one in |30j . Of 
course, such procedure cannot account for the complex- 
ity of non-equilibrium situations encountered in experi- 
ments or simulations [SJ [TTl [TTHTOl [50] . The two-phase 
regime exists only for temperatures T < Tc, and for an 
average density between the value at the liquid and gas 
binodals, as seen in Fig. [2]a. Thus, a piece of fluid ex- 
panding near thermodynamic equilibrium can enter the 
two-phase regime, which amounts to partially undergo- 
ing the liquid-gas first-order phase transition, in only two 
ways: by crossing the liquid binodal, which we call the 
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subcritical case, or by crossing the gas binodal, which we 
call the supercritical case. 

In the subcritical case, which corresponds to the calcu- 
lation of Fig. [l] when it is mapped onto the corresponding 
phase diagram, an expanding monophasic liquid reaches 
the liquid binodal and becomes overstretched. The equi- 
librium configuration for a piece of fluid in the two-phase 
region is a mixture of liquid droplets (of yet undetermined 
size) and gas whose densities are at the binodals for the 
same temperature. We call fragmentation the transfor- 
mation from the monophasic liquid to the non-connected 
cloud of droplets. In the supercritical case, achieved if 
the material is initially heated to higher temperatures, 
a monophasic gas becomes supersaturated after crossing 
the gas binodal, and we call nucleation the process by 
which a certain distribution of liquid droplets is created. 

Figure [2] represents the two cases and the various ex- 
perimental situations that they involve. In Fig. |2]a, we 
show the Van der Waals phase diagram for Al that we 
use in the following, and a schematical representation of 
the sub- and supercritical cases of two-phase expansion 
(arrow 1 and 2). In Fig. [2jb, reproduced from [17 , one 
sees the particle distribution in a 2D MD simulation of 
laser ablation. Due to the inhomogeneous energy depo- 
sition, different types of thermodynamic evolutions are 
seen at a same time, and we use this picture to illustrate 
the various situations of our classification, although not 
following exactly the terminology of the original work 
[31j . In zone I, the dense material is still continuous. In 
zone II, the expanding liquid has undergone cavitation 
of gas bubbles. In the upper zone II and in zone III, 
the liquid is fragmented and the material has entered the 
two-phase regime in the subcritical case. In zone IV, the 
fully atomized, expanding material is likely to reach the 
gas binodal in a later stage, thus corresponding to the 
supercritical case. In the upper zone III, small clusters 
are present, which could originate from fragmentation at 
high temperature (subcritical case) or recent nucleation 
(supercritical case). This is why in Fig. [2] a., where the 
different zones are placed qualitatively, two positions are 
proposed for zone III. 



B. Initial droplet size 

Both cases lead to droplets formation. In order to ini- 
tialize the kinetic model that we present further, it is 
necessary to know the initial droplet size at the onset of 
the two-phase regime. 

In the subcritical case, the overstretched liquid starts 
cavitating (see Fig.[2]b, zone II), which we call the hubbies 
regime and then the bubbles percolate until the liquid 
phase is not continuous anymore (see Fig. [2]b, zone III), 
which we call the droplets regime. We assume that the 
droplets regime starts when the gas and liquid volumes 
are equal: Vg = Vi, which is justified by an argument of 
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FIG. 2. (Color online) (a) Phase diagram of the van der 
Waals EoS for Al (see Sec. Ill) showing the liquid and gas 
binodals (solid lines), and a schematic representation of the 
subcritical (arrow 1) and supercritical (arrow 2) cases of two- 
phase expansion, (b) 2D MD simulation of laser ablation 
with inhomogeneous initial temperature, from [T7], showing 
material in various situations of two-phase expansion, which 
we also locate qualitatively on the schematic classification of 
Fig. [2]a (roman numbers). 



surface energy minimization. 

The mean droplet 's size in a fragmentation scenario 
can be obtained by considering a balance between the 
disruptive inertial forces and the restoring surface tension 
|14) . The model proposed initially by Grady has 
been abundantly validated by MD calculations [HI [171 
[20] in 2D and 3D, and is in very good agreement with 
measurements on He jets [33 . We note that the scaling of 
the mean radius R of the droplet can be obtained by just 
setting to unity the Weber number We = pRv^ /a [34j . 
where a is the surface tension p the liquid mass density, 
and V = rjR the typical velocity difference across a piece 
of fluid of size R. We is the ratio of the surface energy to 
the inertial energy. In any dimension this criterion yields 



We 



R 



PT 



(1) 



Several values of order 1 have been proposed for the pref- 
actor in this scaling law, either from analytical estimates 
(prefactor 15^/^ = 2.47 in [25]), or from fits to MD sim- 
ulation results. In [TS], it was shown that both 3D MD 
results with a homogeneous strain rate 77 and data from 
helium free jets experiments from [33] could be fitted to 
Eq. [1] with the same prefactor, thus validating this law 
over almost 8 orders of magnitude in cluster mass (the 
experimental fragments cover larger sizes than the nu- 
merical ones). 

Concerning the size distribution of droplets resulting 
from fragmentation, MD simulations have shown clearly 
that it is essentially exponential [T4Hl6[[20] . which is con- 
sistent with simple models of entropy maximization (14] . 

By contrast, it is not so clear how to describe the initial 
situation in the supercritical case. This task requires one 
to choose a model for nucleation, or to input results from 
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MD calculations. Nucleation of clusters from a supersat- 
urated vapor is the situation of nucleation whose kinetics 
is the easiest to model theoretically [35] , but still choices 
have to be made [T3] , that are beyond the scope of this 
paper. Any model for nucleation will depend crucially 
on surface tension, so we make the remark here that es- 
timating the surface tension for small droplets is delicate 
because of its enhancement at small sizes [36l EZ] • 



III. MODEL 
A. Van der Waals fluid model 

The kinetic model that we present in Sec. III.C is 
applicable to any EoS. In this paper, for a qualitative in- 
vestigation of the two-phase expansion regimes, including 
kinetic and surface effects, and with emphasis on the late 
time and low temperature limits, we use a van der Waals 
(VdW) fluid model, which for convenience we describe 
first. It is important to note, however, that the VdW 
model is not intended to provide a highly accurate de- 
scription of a fluid, especially in WDM or supercritical 
conditions where ionization and radiation effects can be 
strong. 

With only two parameters, the VdW EoS is the sim- 
plest EoS describing the coexistence of a liquid and a gas 
phase, and has already been used for theoretical studies 
of dynamic two-phase processes [3Sl [32] • All the thermo- 
dynamic functions can be derived from the expression for 
the mean-field potential energy per particle in such fluid: 
U = -\-oo if n > 1/6 and U = —an if n < 1/6, where 
n is the particle density, 6 stands for the incompressible 
volume of the particles, and a represents the mean-field 
attractive energy between them. 

The bulk VdW energy of N particles at temperature 
T is E = N{cyT — an). It can be shown that the spe- 
cific heat Ct, is independent of n and can only depend on 
T [30], so that one has to choose necessarily Cu = |fcs, 
where ks is the Boltzmann constant, if one wants the 
EoS to match the perfect monoatomic gas in the dilute 
limit. Writing the partition function, one obtains the 
other thermodynamic functions. In particular, the pres- 
sure is P = ksT/ (w — 6) — a/v^ where v = 1/n is the vol- 
ume per particle. This expression implies that the isobars 
(isotherms) are a cubic relationship between T and v {P 
and v). Hence, below a certain critical temperature Tc, 
an unstable zone of negative compressibility appears in 
the phase diagram, limited by the two spinodals. We ob- 
tain the equilibrium density of the two stable phases that 
can coexist at certain (P, T) by numerically perform- 
ing the Maxwell construction, which consists of solving 
Pi = Pg (i) and /j,; = /ig (ii) simultaneously, where fi 
denotes the chemical potential and the subscripts / and 
g stand for liquid and gas, respectively, all throughout 
this paper, (ii) is equivalent to vdP = and thus 



If P{v)dV = Pl,g{Vg - Vl) m- 

Introducing the reduced temperature — ksT/lo, 
where Iq — a/b is the latent heat at T = 0, and two 
dimensionless parameters that are small in the low tem- 
perature limit: Vg = b/6 and vi — 6(1 -I- 7), equations (i) 
and (ii) become: 
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It is worth remarking that Tc — 8a/276, so 6'c = 8/27 ~ 
0.3, and therefore one expects that calculations in the 
"low temperature limit" (0^1) should be a good ap- 
proximation as soon as one is not considering the vicinity 
of the critical point. 

Figure [3] gathers the thermodynamic functions of our 
VdW model for aluminum. Figure [3ja shows the nu- 
merical result of the dimensionless Maxwell construction 
where the VdW parameters a = 9.1 x 10~'^^erg.cm'^ and 
6 = 1.85 X 10^^'^cm'^, giving Iq — 3.07ey, have been 
adjusted to fit this material. For that, we impose that 
the VdW liquid density matches the aluminum liquid 
density ri;(r„i) = 5.26 x lO^^cm"'^ at the melting point 
Tm = 933. 5K (= Omeio) [H] and that the VdW latent 
heat (shown in Fig. [sjb) 
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coincides with the experimental value l{Ti,) — 4.88 x 
lO^^^crg/atom for aluminum at the boiling temperature 
Tb = 2792K (= 0.078^0) [31] ■ Note that the critical pa- 
rameters obtained in this way are consistent with the best 
estimates to date, although not very precisely [42] . 

In Fig. [3]a are also displayed the simple, useful ap- 
proximations for ni and ng at lowest orders in 6 that one 
obtains directly from Eq. [2] and [3] : 

mc^lii-e-e^), (5) 
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Note that Eq. [6]is, at lowest order in 9, equivalent to the 
Clausius-Clapeyron formula applied to a perfect gas. We 
also show in figjSjc the approximations at first order in 9 
for the liquid and gas bulk energies per particle 

eg ~ ^fesT , ei ~ ^keT - ^. (7) 

For the surface tension. Van der Waals himself had al- 
ready proposed to model it using density gradients [55] . 
but we have chosen to use a simple formula that is uni- 
versally verified in simple fluids [43 
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FIG. 3. (Color online) Van der Waals thermodynamic func- 
tions for Al, in VdW units. The filled circles represent exper- 
imental data and the empty cirlce the critical point obtained 
from the fits, (a) Liquid and gas densities (solid lines) with 
first (dashed lines) and second order (dotted lines) low T ap- 
proximations, (b) Latent heat (solid lines) decomposed in the 
first (dashed lines) and second (dotted lines) term of Eq. [4] 
(c) Bulk energies per particle, with first order low T approx- 
imations Eq. [t] (dashed lines), (d) Surface tension (Eq. [8|. 



To model aluminum, we fit this formula to the exper- 
imental value a{Tyn) = 1050erg/cm^ [33], as shown in 
fig. [Sjd. Note that in the following, the total liquid en- 
ergy in the cell is Ei = Niei+aSi, where Si is the surface 
area of the droplet. 



B. Kinetic equations. Validity condition 

Our goal is to compute the evolution of droplets in 
cells. In this paper we limit ourselves to the case of one 
droplet in one Lagrangian cell undergoing adiabatic ex- 
pansion. It is a closed system out of equlibrium, and its 
complete description requires the determination of the 
four variables Ni,ni,Ti,Tg. To compute their evolution, 
we need four rate equations: a liquid-gas particle ex- 
change rate, an energy exchange rate, a total energy loss 
rate (work to the outside), and an internal equilibrium 
condition to determine the liquid density. As shown in 
fig. |4ja, the particle fluxes between liquid and gas are di- 
vided between evaporating, condensing, and condensing 
but not-sticking particles. 

The volume expansion V{t) shall later be prescribed by 
a global hydrodynamic code. For our study, we assume 
cylindrical symmetry and we use a simple model behavior 

[13 [H] 
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where Vq = Lq is the inital cell volume and rjz = 
{dvz/dz)t=o and 77,- = (dvr/dr)t=Q are the axial and ra- 



dial strain rates, respectively defined as the initial veloc- 
ity gradients in the beam direction and the target plane. 

The evaporation and condensation fluxes are computed 
using the standard Hertz-Knudsen formulas 
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where m is the particle mass and n* (Ti , R) is the equi- 
librium gas density for a droplet at temperature T; and 
of radius R. To estimate n*, we use the Kelvin equa- 
tion, which describes the increase of the equilibrium va- 
por pressure surrounding a droplet due to surface tension 
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The Kelvin equation is approximate because its deriva- 
tion assumes a perfect gas. Also, we use a constant value 
for CT, thus neglecting its increase at small radii [36 l 137 ) . 
Still, this approach is probably not too bad after all [3H] , 
and satisfactory enough for our qualitative purpose. 
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FIG. 4. (Color online) The "droplet-in-cell" kinetic model, 
a) Sketch of the kinetic fluxes, b) Sticking coefficient P cal- 
culated with formula from [49] and our aluminum VdW pa- 
rameters. 



Considering mass conservation, and combining the two 
fluxes of Eq. |10[ the particle exchange rate equations are 



d{Ni+Ng) 
dt 



= 



dNi 
dt 



= /?(-*vap + $cond)5;, (12) 



where Si is the surface area of the droplet and < /3 < 1 
the sticking coefficient that is usually assumed of order 
0.5. A recent study [33 has proposed a simple expression 
for /3 that is in good agreement with MD calculations for 
several simple fluids. This expression depends only on 
the ratio of the molecular volumes in the liquid and vapor 
phase: /3 = (1 - (wi/wg)^/^) exp(-l- 



— — t) that we 

plot for our VdW model for Al in Fig.|4]b. Note however 
that, for the calculations presented in this paper, we use 
a constant value (3 = 0.5. 

Concerning the energy fluxes, the first equation comes 
from the assumption of adiabatic expansion of the cell : 



djEi + Eg 
dt 
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In this global energy loss rate we have neglected three 
terms that could be added in the near future. The first 
one is heat conduction between cells. This term may 
play a role, but it cannot be very important as we are 
considering a supersonic flow (see Sec. I). The second 
neglected term is radiation. Radiation becomes indeed 
the dominant cooling mechanism at long times, as we see 
later, but it is negligible for the initial dynamics, so the 
approximation is reasonable, because our purpose in this 
paper is to study the expansion in a time range where the 
two phases are interacting and the system is not just a 
collection of isolated clusters flying in vacuum. The third 
neglected term is thermionic emission. One expects elec- 
trons to be thermally emitted from the droplet, taking 
away some energy. Non-neutral effects are totally absent 
from our model, but we expect that the associated cool- 
ing rates will be small compared to the adiabatic and 
evaporative coohng rates [50] . 

The energy exchange rate between the liquid and the 
gas has contributions from the three fluxes of Fig. |4]a. 
The contribution of the colliding but non-sticking parti- 
cles can be described with a flux proportional to the tem- 
perature difference Tg — Ti, with a relaxation coeficient 

< a < 1 (see e.g. [5T] for more discussion). For the 
condensing gas particles, we make the simplest assump- 
tion, that each of them brings into the liquid the average 
gas energy per particle Cg. For the evaporating particles, 
we assume that the energy they individually take away 
from the liquid depends only on the liquid state. We note 
it e* and define it as the energy of a virtual gas particle 
that would be in equilibrium with the droplet of radius R 
at temperature T;. For the VdW fluid, those energies are 
Cg = c^Tg — aug aud e* — CyTi — an*g{Ti, R), respectively. 
Note that our definition of e* is totally analogous to the 
Hertz-Knudsen derivation of the mass evaporation rate 
of Eq. [12] We finally get the exchange rate 

= [/5(-e;«'vap + eg$co„d) 

+ (1 - P)ac,{Tg - TO^cond] Si. (14) 

Our set of kinetic equations is fully consistent in the 
sense that at equilibrium both mass and energy fluxes 
between the droplet and the gas are in balance. In partic- 
ular, the average energy that a liquid particle takes from 
the rest of the liquid to evaporate is e* — e;. This term, 
which for the VdW fluid is — a(n* — ni), corresponds ex- 
actly to the latent heat per particle for an arbitrary fluid 

1 — {cg — ei) — Pi,g{vg — vi), without the second (work) 
term, which is expected since the latent heat is an en- 
thalpy and we are here dealing with energy exchanges at 
constant volume. 

To our knowledge, our set of rate equations is an orig- 
inal model for the exchanges between a droplet and its 
vapor. Other sets of kinetic equations can be found for 
analogous systems (see e.g. 1S51[S3]), but they do not cor- 
respond to the purely kinetic regime that we are consid- 



ering, because they deal with larger droplets (i? > 1/xm) 
and longer time scales, more relevant to the fields of com- 
bustion or atmospheric sciences, so they need to combine 
the kinetic approach with the more classic hydrodynamic 
theories of droplet evaporation [S3] . 

Our model for WDM situations is simpler, because we 
do not distinguish in the gas a Knudsen layer vs a hy- 
drodynamic layer. We assume that our computing cells 
are small enough that the gas density inside them is con- 
stant. The variations over the whole flow shall instead 
be treated by the global hydrodynamic code that deter- 
mines the expansion of each kinetic cell. The kinetic and 
phase change processes in our description are driven by 
the hydrodynamic expansion, therefore, the validity con- 
dition of our model is that the initial cell size should be 
much smaller than the initial dimensions of the expand- 
ing material : 

Lo^Sr,6z, (15) 
where Sz is the initial foil thickness and Sr the heating 



beam diameter. If Eq. 15 is verified, the global hydrody- 
namic treatment is correct, with gradients properly re- 
solved. Eq. [15] is verified in the standard situations we 
consider, but can break down if the initial droplet or cell 
size is not small enough compared to the sample size. 



C. Equilibrium condition between droplet and gas 

In order to get a closed system of equations for particle 
and energy fluxes, we still need one assumption. Our 
model is a priori out of thermal equilibrimn (T; ^ Tg), so 
the density of the liquid is not determined yet. It seems 
reasonable to assume pressure equilibrium between the 
droplet and the gas [T3], because we expect that a few 
collisions are sufficient for the droplet to "experience" 
the gas pressure, and adjusting the liquid pressure to it 
requires only a small density change, due to the very low 
compressibility of the liquid. 

Due to the droplet curvature, the pressure equilibrium 
condition is the Laplace equation 



P,-Pr, 



2a 
li' 



(16) 



An exact numerical implementation of Eq. 16 is difficult, 
because it requires to solve a non-linear system at each 
timestep in order to determine the liquid density given a 
certain set of values {V, Ni, Ng, Ei, Eg}. 

To simplify the condition, one can approximate the 
liquid density by the equilibrium value n*{oo) for a flat 
interface (i? — >■ oo), but this is wrong for two reasons : 
first, because due to the fast expansion, the gas pressure 
is lower than the saturation value corresponding to the 
liquid temperature, and second, because of the Laplace 



compression term of Eq. 16 Running our model, we have 
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checked that this raw approximation leads to important 
inaccuracies in the calculation of the pressure, especially 
at low temperatures where the Laplace term becomes 
dominant. Still, these errors do not cause important dis- 
crepancies in the global description of the droplet evolu- 
tion, which we attribute again to the low compressibility 
of the liquid. 

For more accuracy, we choose to compute the liquid 
density in perturbation from the flat equilibrium value : 

n. = nrM(l + ^), (17) 

where Ki{Ti) — ni{dP/dni)Ti is the isothermal bulk 
modulus of the liquid that we compute directly from the 
VdW EoS, and AP = 2a /R- (P, (n* (oo), T,) - Pg) is the 
pressure correction that we compute using Eq. |16| As we 
show in the next section, this perturbative approach of 
the pressure equilibrium condition is very satisfactory. 

With Eq. [9p7l we have a complete kinetic model for 
the evolution of a droplet and its vapor in a cell. In the 
following, we combine it with the VdW EoS to study the 
different regimes of two-phase expansion. 



IV. RESULTS 
A. NDCX II reference case (subcritical case) 

The reference case envisioned as an upcoming experi- 
ment on the NDCX II machine at LBNL consists of heat- 
ing an aluminum foil of thickness Sz — 3.5/im with a short 
pulse of ions of duration Cheating S Ins, which makes the 
picture of rapid heating roughly correct [7]. The beam 
profile is taken as a uniform disk of diameter 6r = 1mm. 
Initial temperatures up to leV are predicted for the ex- 
pected beam fluences [24] . 

In Fig.[5j we present the numerical output of the model 
for a cell containing a droplet and gas initially at equi- 
librium at To = 8000K with Vi ^ Vg ^ V/2, because 
this corresponds to the onset of the "droplets regime" 
(see Sec. II. B). As we mentioned previously, we make 
the crude assumption of a flow with linear speed profile 
and outward expanding speed vq = 3cq ~ 5.0kni/s on 
both sides z > and z < 0, where the initial sound 
speed Cq is estimated roughly as the thermal velocity 
Vth{To) = -^/fesTo/m. Then, the strain rates in Eq. 9 
are simply r]z = 6cq/Sz and rjr = 6cg/5r. We dis- 
play a full 3D case (solid lines) and a ID case (dashed 
lines) where r]r — 0. Here the hydrodynamic time is 
ihydro = ^/Vz — 0.37ns. The time t^o ~ — 107ns 
can be considered as the time of the crossover to the 
3D regime of expansion. The calculation is carried out 
with a ^ (3 — 0.5 and we use the variable u = ln{t) to 
span a wide temporal range. The result is displayed up 
to t = 100/is because at this time the front has traveled 



over about 50cm, which is comparable with the size of 
an experiment. The initial droplet radius Rq = 25.4nm 
is estimated using Eq. [l] and is consistent with observa- 
tions for similar initial temperatures [13] . This size is the 
mean size of the liquid fragments so we are considering 
the evolution of a typical droplet. 




t (ns) n (cm ^) 

FIG. 5. (Color online) Droplet and gas evolution in the NDCX 
II reference case. Initially the droplet of radius -Ro ~ 25.4nm 
and the gas have equal volumes and are in equilibrium at 
To = 8000K. All liquiwd (gas) quantities are labeled with I 
(g) . Plotted are the time evolution of (a) the particle numbers 
and (b) the temperatures, for a ID (dashed lines) or a full 
3D expansion (solid lines), (c) Pressure evolution in the 3D 
case. The pressure difference computed (dashed lines) and 
expected from Eq. |16| (dotted lines) are indistinguishable, (d) 
Trajectories in the phase diagram for the ID (dashed lines) 
and 3D (solid lines) cases 

At early times {t < 10ns), a fraction of the liquid is 
evaporated (Fig.[5]a). However, this process saturates at 
a time tmin, after which the liquid fraction starts growing 
slowly. Then, in the purely ID case (dashed lines), the 
droplet continues to grow steadily. In the more realis- 
tic situation however (solid lines), the droplet evaporates 
again when the 3D regime sets in, at times t > t^jj. 

In Fig.[5]b, one sees that, almost instantaneously after 
heating {t < lOOps), a temperature difference AT = Ti — 
Tg is established between the gas and the droplet, and 
remains roughly constant throughout the expansion in 
the ID case. On the contrary, in the 3D case, Tg drops 
quickly to almost around tso, whereas T; decreases 
slowly to a value around 1600K. 

In the phase diagram trajectories [Fig.[5]d], we see that 
in both cases the liquid density remains very close to 
the equilibrium value. By contrast, the gas density is 
clearly below the binodal in the ID case, and in the 3D 
case it dives deep into non-equilibrium (supersaturated) 
conditions. 

In Fig.[5]c, we check the pressure equilibrium condition 
in the 3D expansion case. One cannot distinguish the 
pressure difference in the computed evolution (dashed 
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lines) that uses Eq. 17 from the theoretical value of Eq. 16 
(dotted lines), as the agreement is better than 2% over 
the whole simulation range. The increase of Pi at long 



times is due to the Laplace term (Eq. 16). 

Clearly, from the NDCX II example, two different 
regimes can be identified. The first one, where the tem- 
perature difference is small, and remains constant, is 
a quasi-thermalized regime. In this regime the droplet 
grows. The second one, where the gas becomes much 
colder than the drop, is a non-thermalized regime. In 
this regime the droplet evaporates again, as if it were in 
vacuum. We now discuss the various regimes. 



B. Thermalization condition, quasi-thermalized 
regime 



the gas. With the approximation i?o — -R and neglecting 
the prefactors of order unity, we obtain the simple scaling 
that is expected to be valid for any EoS : 
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Let us find a thermalization condition. In our equa- 
tions, the energy is extracted from the system only by 
the adiabatic expansion of the gas (Eq. \l3^ and the gas 
quenching is then transmitted to the liquid via the liquid- 



gas energy exchange term (Eq. 14). Therefore, we should 
compare those two energy fluxes to find the thermaliza- 
tion condition. 

Let us assume a small temperature difference AT/T ^ 
1, so that we are in the quasi-thermalized regime of ex- 
pansion, as in the ID case of Fig. [5] From Fig. [5]a, one 
sees that Ni and Ng are almost stationary if T; ~ Tg. 
Hence, let us make the approximation ^vap — 'i'cond 
(more precisely, |$vap - $cond| < "J-cond)- 

The ratio between the two energy fluxes can be esti- 
mated as follows. Let us consider a cell containing fixed 
numbers Ni and Ng of liquid and vapor atoms, respec- 
tively, and let x = Ni/Nt^t be the liquid fraction in the 
cell, where A'tot — Ni + Ng. In the low T limit, a small 
energy change can be written dEg = Ng^ksdTg for the 
gas and dEi = Ng^ksdTi for the hquid, according to 
Eq. [T] Noting dEtot = dEi -\- dEg the total energy lost by 
the cell, we define ^ = dEi/dE^ot- Requiring stationary 
AT (i.e. dTi = dTg), we obtain ^ = 5x/{2x + 3). 

In the low T limit, the adiabatic cooling of the gas 
implies: dEtot/dt = —PgdV/dt — —rigkBTgfjVo, where 
we define the volume strain rate 77 = d{V/Vo) /dt. Note 
that, in the ID expansion regime, fj = rj^, whereas in 
the 3D expansion regime, for t ^ t^jj, fj :^ Sri^rj^t^ and 
diverges. The power transferred from the liquid to the 
gas is: dEi/dt = ngy/kBTg/2'KmSxCvAT, where x = 
l3 + {1 — l3)a. To get this expression we have computed 
the contributions from the three terms in Eq. [Ml and 



used $ 



vap 



^cond- Expressing S = AnR^ and Vq 



2 X |7ri?Q, the balance between the fluxes dEi = ^ dEt. 



finally yields 



AT 4V27r 



7?3 



X 9 vthiTg)R^ 



(18) 
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(Color online) Test 
Tsl (a) Time 



of the thermalization for- 
evolution of the ratio A = 



co/ (AT/T) of the temperature difference computed 
|18| to the numerical model in the NDCXII refer- 
ence case, in purely ID (dashed line) and full 3D (solid line) 
expansion, (b) A in the ID expansion at ti = Ins (solid sym- 
bols) and at t2 = lOps (open symbols) with p —0, 0.2, 0.4, 
1 (circles), a —Q, 0.2, 0.4, 1 (up triangles), Sz =1, 2, 4, 8/im 
(squares) and -Ro =5, 10, 20, 40nm (down triangles). 



In Figure [6] we check the validity of Eq. [18] comput- 
ing the ratio A = (Ar/r)thoo/(Ar/r) of the theoretical 
temperature difference (Eq. |l8| ) to the result of the full 
numerical calculation. To evaluate Eq. 18 we take the 



values of fj, R, Tg and x from the result of the numeri- 
cal simulation. In Fig. [6]a, we show the evolution of the 
ratio A in the NDCXII reference case (same calculation 
as Fig. [sj). In the full 3D expansion case (solid line), 
the prediction becomes bad (error larger than 100%) at 
t ~ t^D, which is expected since the volume expansion 
rate fj diverges in 3D. In the purely ID expansion, after 



18 is obtained 



the first ns, one sees that Eq. 18 is accurate within 20%. 
More precisely, the error has the same sign as the deriva- 
tive dNi/dt and vanishes when the droplet is stationary, 
at t = tinin, which is expected since Eq, 
with the assumption dNi/dt — 0. 

In Fig. [6|b, we show the ratio A at t 

those that are relevant to Eq 



1 — Ins and at 
Ifis for the same parameters, but varying one by one 

/3, a, 5z (in order to 



18 



where we note Wth(Tg) — ^JkBTgjm the thermal speed in 



vary 77) and Rq. The analytic formula overestimates (un- 
derestimates) AT/T in all cases at ti (^2), and for both 
cases the error is larger when the expected (AT/T)theo 
is larger, which is natural since Eq. [18] holds in the limit 
of small AT/T. Interestingly, the point /? = is sepa- 
rate from the others at both times, and is closer to 1, 
which is not surprizing since 13 — Q means no particle ex- 
change, and this again confirms that the main source of 
inaccuracy of Eq. 18 is a non-zero value of dNi/dt. The 
prediction could thus be refined if this effect was taken 
into account, for example using the analytical results of 
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the next section. At early times, one can see also that 
the error is larger for droplets of radii smaller than 5nm. 
This effect is switched off if we set ct = 0, confirming that 
it is caused by surface effects. 

It is also possible to rewrite Eq. [19] in terms of the ini- 
tial conditions: considering r] ^ Cq/6z, one gets the very 
simple scaling law: AT/T ~ Rq/6z This last expression 
is only valid in the case of a ID linear expansion where 77 
is constant. In this case, one sees using Eq. [ijthat AT/T 
is expected to decrease slowly when the initial sample 



size increases 



AT 



cx 773 (X Sz 



(20) 



For larger samples, the thermalization will be better even 
though the droplets are bigger. This justifies again that 
in the limit of large samples and slow expansions, an 
equilibrium hydrodynamic description becomes valid. 
In summary, Eq. 18 is expected to be always a good 



estimate in the quasi-thermalized regime, and Eq. |19| can 
be considered as a universal criterion to delimit the quasi- 
thermalized regime. 



in fig. [5] Note that, at long times, and independently 
from the EoS, droplets will always grow in a thermalized 
situation. This is due to the fact that adiabatic expansion 
of a perfect gas is an algebraic trajectory in phase-space 
(T cx p^/'^), whereas Clausius-Clapeyron law predicts an 
exponential curve for the gas binodal, meaning that the 
gas in a two-phase expanding cell will always tend to 
saturate and make the liquid fraction grow. 
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C. Fully thermalized regime 

In the previous section we could distinguish the 
regimes of quasi-thermalized versus non-thermalized ex- 
pansion. From Eq. 18 it is clear that the quasi- 
thermalized regime will become quickly invalid after t^o , 
because the volume strain rate 77 diverges. Nonetheless, 
in the early times of expansion, or if one is interested in 
systems of large radial extent, it is worth studying the 
limiting case of a fully thermalized flow. 

In this perspective, let us assume Ti — Tg ^ T. Again, 
we look at the low T regime, which becomes valid very 
early in the expansion process. Using the first order ap- 
proximations in the VdW model ri; ~ (l — 9)/b and Ug ~ 
(Eq. [5] and [6]) and neglecting the surface energy term, 
we write the total energy Etot = Ntot{cvT — x{l — 9)a/b), 
where x = Ni/Ntot is still the liquid fraction in the cell. 
The total energy change dEtot — ~PgdV becomes, at 
first order in 9: -(1 - x) 9 = + x) d0 - dx. Noting 
that dV ~ dVg, we convert 9 d{\n{V)) ^ d9 - d9/9 using 
the low T approximation Eq. |6] and find finally 



dx = 



1-x 



d9. 



(21) 



It is easy to push the approximation to higher orders in 9, 
but Eq. [21] already allows one to get good insight into the 
evolution of the droplet. One sees that the droplet will 
be stationary at a temperature satisfying 6'min = §(1 ^ 
x) corresponding to the time imin already mentioned, it 
will evaporate before this point, for temperatures 9 > 
Omim and grow after it, for 9 < ^min- This sequence is 
in agreement with the NDCX II reference case shown 



FIG. 7. (Color online) Thermalized evolution for a droplet 
of initial radius Ro = 25nm and initial temperatures To = 
7000 K to lOOOOK (solid lines). Temperature (a) and liquid 
fraction (b) evolution versus the cell expansion, (c) Full nu- 
merical calculation compared to solution of Eq. 
lines) started at V/Vo = 10. 



21 



(dashed 



In Figure [7j we show the exact numerical calculation 
of the thermalized evolution of a droplet whose initial ra- 
dius is Rq — 25nm and for initial temperatures ranging 
from 7000 to lOOOOK. In the thermalized case, the time 
evolution is irrelevant, that is why in fig. [Tja and b the 
temperature and liquid fraction are plotted as a func- 
tions of the volume expansion ratio V/Vq- An expansion 
remaining in the thermalized regime over 10 orders of 
magnitude volume expansion is unrealistic in the case of 
NDCX II, but for generality it is interesting to study this 
limiting trend. 

In Fig.[7]c, the numerical liquid fraction versus temper- 
ature is compared to the solution of Eq. |21[ starting when 
the volume has expanded by one order of magnitude. 
The analytic approximation is accurate within 20%. This 
shows that Eq.[2l]can be used to make good estimates of 
the asymptotic growth of the liquid fraction in the ther- 
malized case, which can be, e.g. for Tq = 9000K, of a bit 
more than 10%. 

This growth of the liquid fraction in the thermalized 
regime is a rigorous upper bound for droplet growth. In 
the opposite regime, one can get the reciprocal upper 
bound for droplet evaporation. 
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D. Non-thermalized regime: evaporation in vacuum 

If the gas expands too fast for thermalization to oc- 
cur, one expects the droplet to evaporate as if it were in 
vacuum. The corresponding limit consists of assuming 
^cond — 0. Obviously, in this case the droplet can only 
lose particles. Moreover, the evaporation of the droplet 
is maximal in this regime, because, if there was thermal- 
ization via collisions with a gas, colder than the drop, 
it would be a way for the droplet to lose energy with- 
out losing particles. Also, because the vapor pressure 
decreases exponentially with temperature, one expects 
the evaporation in vacuum to slow down fast. However, 
independently from the kinetics, it is clear that there 
must be some upper bound to the evaporation of a drop. 
Indeed, as every evaporating particle takes away energy 
from the droplet (the latent heat), the droplet gets colder 
and colder, until the evaporation is "frozen" , a strict limit 
being that T cannot become negative. 

Let us find analytic expressions for the maximal evapo- 
ration of a droplet whose energy is noted £"/ . Considering 
the evaporating particles and using Eq. [T4j the energy 
loss can be written dEi = e*gdNi. On the other hand, 
considering the liquid, and neglecting the surface energy 
term, one has dEi ~ eidNi + Ni{dei/dTi)dTi. Equating 
those two expressions, one finds, for the VdW fluid : 



dNi 



(c„ - a{dni/dTi)) 



dTi. 



a[ni 



(22) 



Using the development of ni at first order in 9 (Eq. [5]) , 
one can integrate Eq. [22] from initial Tq and Niq, yielding 
at the final T/ the remaining fraction 



exp 



4 ^ 



9f) 



(23) 
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FIG. 8. (Color online) Evaporation in vacuum. Time evolu- 
tion of the temperature (a) and the non-evaporated fraction 
(b) for a droplet of initial radius _Ro ~ 25nm and initial tem- 
peratures To = 3000 to 9000K (solid lines). The dashed lines 
on (b) are the numerical integration of Eq. 
lines correspond to Eq. [23] 



22 



and the dotted 



Figure ]8] shows the exact numerical result for the evap- 
oration in vacuum of a droplet whose initial radius is 
25nm, and initial temperatures ranging from 3000 to 



9000K. The volume expansion is not relevant here, so 
the variables are displayed as a function of t only. In 
Fig. l8]a and [8]b, one sees that for all initial tempera- 
tures, the number evaporation and cooling curves of the 
droplet follow a same asymptotic behavior, which is in- 
creasingly slow at long times. In Fig. [8]b, the numerical 
integration of Eq. 22 is shown for each Tq (dashed lines) , 
with the final T; taken from the full numerical solution. 
The agreement with the final evaporation ratio is excel- 
lent, showing that surface effects play a negligible role. 
We have checked that surface effects cause an overestima- 
tion of the maximal evaporation of less than 10% even for 



droplets of initial radius Inm. The approximated Eq. 23 
is also displayed for each Tq (dotted lines), taking here 
also the final T; from the full numerical run. One sees 
that it predicts the good limit for evaporation within 5% 
for initial temperatures up to 7000K. This is very satis- 
factory because the non-thermalized regime is expected 
to be valid only in the late times of expansion so the first 
order low T approximations should be very accurate. 

Let us now discuss the onset of the radiative cooling 
regime. At long times and low temperatures, the par- 
ticle evaporation and the evaporative cooling rates de- 
crease exponentially (Eq. [6]) , whereas the radiative cool- 
ing rate is algebraic (oc T^). Therefore thermal radia- 
tion becomes the dominant cooling mechanism at long 
times. Within our model it is not difficult to express the 
temperature Trad below which radiative cooling becomes 
dominant over evaporative cooling [55,. Using our VdW 
parameters for Al, and assuming an emissivity e — 0.2, 
we find Trad — 1740K. For Si nanoparticles, the same 
estimate yields Trad — 2190K. This is consistent with 
the measurements reported in [T^] where the cooling of 
Si nanoparticles formed by laser ablation is found to be 
well explained by radiation at expansion times 5 — 150/zs 
and for temperatures below 2000K, although the evapo- 
rative cooling rates that we obtain within our model are 
significantly larger than the estimates they report. As 
an example, for Si at 2000K, our crude model predicts a 
radiative cooling rate of ~28K//is, while the evaporative 
coling rate in vacuum is still ~6.5K//is. Note however 
that the rate we compute is a strict upper bound. 



E. Supercritical case: nucleation and growth of 
liquid droplets 

The last case to consider is the supercritical case, where 
the material expands first as a supercritical fluid, and 
enters the two-phase region of the phase diagram crossing 
the gas binodal, as a supersaturated gas. In this case 
nucleation of droplets may occur. 

We do not propose a model for nucleation but we note 
that it has been reported [T3] that supersaturation of the 
vapor does not reach high values and that above a cer- 
tain threshold value nucleation is very sudden, due to 
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the exponential dependence of the nucleation frequency 
on the supersaturation ratio |35j . Then, our model is 
expected to describe correctly the subsequent evolution 
of the clusters. In particular, we expect that thermaliza- 
tion will depend on the droplet radius and the volume 



strain rate and that Eq. 18 will be a good estimate in 
the quasi-thermalized regime. If clusters and gas are in 
thermal equilibrium, we expect Eq.[2l]to be valid as well. 

Note that nucleation may also happen in a subcritical 
expansion scenario, if the gas becomes very supersatu- 
rated at long times, as can be seen for the reference case 
in Fig.[5]d. In this situation, close to the non-thermalized 
limit, condensation on the existing droplets is too slow, 
and nucleation of new droplets may happen even with 
liquid clusters already present in the plume. 



V. CONCLUSION 

In conclusion, we have studied droplet evolution and 
thermalization conditions with an original, simple kinetic 
model based on a consistent set of rate equations for mass 
and energy exchanges. Using the vdW EoS as a test-bed, 
we have demonstrated that such a kinetic treatment is 
able to bridge the gap between the molecular and equi- 
librium hydrodynamic approaches that have mainly been 
used so far. Most of the results of this work are general 
and should be extendable to any EoS with which the ki- 
netic equations are used. 

In particular, the main output of our study is to iden- 
tify the different regimes of two-phase expansion. On one 
side, the quasi-thermalized case and its limit, the fuUy- 
thermalized case, on the other side, the non-thermalized 
case. To distinguish the two situations, we identify a lo- 



cal thermalization condition (Eq. 18 ) which depends on 



the droplet radius R, the volume expansion rate rj, the 
gas temperature Tg, the liquid fraction a;, and the ki- 
netic parameters a and /3, but it can also be traced back 
to the initial conditions: sample thickness Sz and initial 



temperature Tq (Eq. 20). Eq. 19 is a simpler alternative 
to Eq. 18 that involves only the initial droplet radius i?o, 



fj and Tg and gives the scaling expected for any EoS. 

Due to the crossover from ID to 3D expansion at the 
time tsD, the expansion is expected to take place in the 
quasi-thermalized regime at early times, at least in the 
NDCXII reference case and for similar parameters, but 
at long times the non-thermalized regime is almost in- 
evitable. Eq. [18] shows that this is only a dimensional 
effect, driven by the divergence of fj in 3D. 

In the quasi-thermalized case, our study shows that 
the relative temperature difference (T; — Tg)/Ti remains 
almost constant throughout the expansion. Eq. [18] is de- 
rived assuming no net particle exchange, so only kinetic 
energy terms (but no latent heat) are involved, which 
makes it suitable for generalization to other EoS. The 
predictions of Eq. |18[ and, more generally, the validity of 



the whole droplet evolution model, could also be tested 
with MD calculations and experimental measurements. 

In the fully thermalized case, droplets can grow (mod- 
erately) at long times and we give an approximate for- 
mula for the VdW fluid (Eq. 21). In the opposite case 



of a fully non-thermalized flow, one can find a strict up- 
per bound for the evaporated fraction at a given final Ti . 



For the VdW fluid, this limit is Eq. [22] In both limiting 
cases, simple implementations of the kinetic model, e.g., 
with the VdW EoS, can be used to make estimates and 
provide upper bounds for the droplets evolution. Note 
that, in all the cases we have studied, the droplets never 
grow or evaporate very much from their initial situation. 

For the moment, the model we have presented is lo- 
cal, but in the future it could become part of a larger 
hydrodynamic code that will treat many lagrangian two- 
phase cells containing droplets and gas. The extension 
of our model to several droplets in one cell can be done 
easily. For more realistic simulations, it will probably 
be necessary to take into account other effects that are 
not two-phase phenomena and that we have left aside, 
such as radiation, which can be non-blackbody in the 
early stages of expansion, but also, thermal conduction 
between cells, and thermionic emission. 

Before our kinetic model can be used to compute 
droplets evolutions at the core of a global comprehen- 
sive code, additional modules are needed to initialize the 
two-phase regime. In the subcritical hydrody- 
namic code and some model for fragmentation is required 
to determine the droplets distribution at each location, 
and possibly their velocities, which may differ from that 
of the expanding gas. In the supercritical model 
for nucleation is needed. In any situation, the thermal- 
ization condition (Eq. 18 or 19) can be used to check 



wether the two-phase computing cell can be treated as 
an equilibrium cell or if a non-equilibrium treatment is 
required. The reason being of course that an equilibrium 
(thermalized) description is much easier to implement. 

Finaly, we have investigated the role of surface effects 
in different cases. Surface tension is expected to play an 
important role for droplets of radii R < R^ = a/ksTni. 
This can be seen from Kelvin equation or considering the 
radius at which the surface energy becomes comparable 
to the kinetic energy per particle in the liquid. With 
our VdW parameters for aluminum, R^, increases from 
0.8nm at T = lOOOOK to 64nm at T = 2000K. Surface 
effects are thus increasingly important in the late stages 
of expansion, at low temperature and for the smallest 
fragments. This is also why a careful treatment of the 
supercritical case of in-flight nucleation is more difficult 
and remains to be done in order to complement this work. 
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